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1      Introduction 

Intersection  problems  are  among  the  fundamental  topics  in  computational  geometry 
and  geometric  modeling.  They  are  also  considered  important  in  computer  animation, 
physical  based  modeling,  robotics,  and  computer  simulated  or  virtual  environments.  In 
most  of  these  applications,  the  objects  are  undergoing  motion  governed  by  dynamic 
constraints  or  external  forces.  It  is  important  to  determine  contact  or  interference 
between  any  pair  of  objects  in  a  dynamic  environment  and  generate  a  collision  response 
to  it.  This  problem  has  received  a  lot  of  attention  over  the  last  few  years  in  all  these 
areas. 

The  problems  of  contact  determination  and  interference  detection  have  been  ex- 
tensively studied  in  different  areas.  The  literature  in  computational  geometry  consists 
of  a  number  of  theoretically  efficient  algorithms  for  polyhedral  objects  in  static  envi- 
ronments [1,  2,  8,  9,  10,  11,  12,  16,  34].  There  are  a  number  of  algorithms  with  good 
asymptotic  bounds,  however  their  practical  utility  is  not  clear  since  many  of  them  are 
not  implemented  in  a  realistic  environment. 

Many  algorithms  are  known  for  intersection  between  curved  objects  represented 
as  algebraic  surfaces  or  piecewise  spline  models  in  geometric  modeling  [13,  24,  31, 
23,  35].  However,  they  are  targeted  towards  robust  computations  of  the  intersection 
curve  between  such  models  for  CSG  operations  and  boundary  computations.  They 
are  usually  specialized  to  static  environments  and  are  relatively  slow  for  the  kind  of 
performance  desired  in  computer  simulated  environments 

Similarly,  the  problem  of  collision  detection  has  been  addressed  in  robotics  liter- 
ature as  well.  However,  the  main  purpose  of  applications  has  been  on  planning  a 
collision-free  trajectory  between  obstacles  [7,  14,  26,  29,  40].  This  is  different  from 
applications  in  physical  based  modeling  and  virtual  environments,  where  the  motion 
is  subject  to  dynamic  constraints  or  external  forces  and  cannot  typically  be  expressed 
as  a  closed  form  function  of  time.  . 

The  simplest  algorithms  for  interference  detection  in  dynamic  environments  are 
based  upon  using  bounding  volumes  and  spatial  decomposition  techniques.  Typical 
examples  of  bounding  volumes  include  cuboids,  spheres,  octrees  etc.  and  they  are 
chosen  due  to  the  simplicity  of  finding  interference  between  two  such  volumes.  These 
bounding  volumes  work  very  well  in  performing  "rejection  tests",  when  two  objects 
are  far  apart.  However,  once  the  two  objects  are  in  the  vicinity  of  each  other,  spa- 
tial decomposition  techniques  based  on  subdivision  are  used  to  solve  the  interference 
problem.  This  can  become  rather  slow  in  practice  as  highlighted  by  Hahn  in  [21]. 

In  computer  animation,  algorithms  for  polyhedral  objects  of  complexity  0(n2m2), 
where  m  is  the  number  of  polyhedra  with  n  vertices  per  polyhedron,  are  described  by 


Moore  and  Wilhelms  in  [33].  For  convex  polytopes,  linear  time  algorithms  based  on 
linear  programming  by  Megiddo  and  Seidel  [32,  39]  and  tracking  closest  points  between 
incremental  motion  by  Gilbert  and  etc.  [20]  have  been  proposed.  However,  they  do  not 
make  use  of  coherence  in  their  methods  between  successive  instances.  More  recently, 
geometric  coherence  has  been  utilized  to  devise  algorithms  based  on  local  features 
for  convex  polytopes  by  Baraff  and  Lin  [3,  27].  This  has  significantly  improved  the 
performance  of  interference  detection  algorithms  in  dynamic  environments.  As  far 
as  curved  models  are  concerned,  algorithms  based  on  interval  arithmetic  have  been 
proposed  by  Duff,  Herzen,  Snyder,  etc.  [15,  22,  17].  They  are  relatively  slow  in 
practice  and  expect  that  the  motion  can  be  expressed  as  a  closed  form  function  of 
time.  However,  no  practical  and  efficient  interference  detection  algorithms  are  known 
for  general  geometric  models  in  dynamic  environments. 

In  this  paper,  we  present  efficient  algorithms  for  contact  determination  and  inter- 
ference detection  between  two  or  more  objects  undergoing  rigid  motion.  The  models 
include  polyhedral  models  and  curved  objects  whose  boundary  can  be  represented 
in  terms  of  algebraic  sets.  This  includes  algebraic  surfaces  and  NURBS  models  fre- 
quently used  in  computer  graphics  and  geometric  modeling.  No  assumption  is  made 
on  the  motion  of  the  objects.  At  each  instance,  their  position  is  described  using  a 
transformation  matrix  with  respect  to  the  origin.  For  convex  polytopes,  the  algorithm 
keeps  track  of  closest  features  based  on  convexity  of  polyhedra  as  described  in  [27] 
and  [28].  Concave  polytopes  are  decomposed  into  convex  polytopes  and  represented 
by  sub-part  hierarchical  tree. 

For  curved  models  we  use  a  combination  of  hierarchical  representations,  algorithms 
for  polytopes,  and  global  and  local  methods  for  solving  systems  of  polynomial  equa- 
tions. In  particular,  we  use  the  control  polytopes  of  the  objects  and  their  sub-patches 
along  with  sub-part  hierarchical  tree  representation  to  describe  them  for  top-level 
collision  detection.  Once  the  control  polytopes  come  into  contact  with  each  other, 
we  algebraically  formulate  the  problem  of  contact  determination  and  interference  de- 
tection. The  problem  is  reduced  to  solving  a  system  of  overconstrained  algebraic 
equations.  We  initially  use  global  algebraic  methods  based  on  resultants  and  matrix 
computations  to  compute  the  solutions  in  the  corresponding  domain  of  interest  and 
update  these  solutions  using  local  numerical  methods  as  the  objects  undergo  motion. 

One  of  the  main  contribution  of  this  paper  is  to  present  a  methodology  which  cap- 
tures the  temporal  and  spatial  coherence  between  successive  checks  to  achieve  real-time 
collision  detection  for  polyhedral  objects  and  to  determine  the  contact  information  be- 
tween general  curved  models  efficiently  using  the  algebraic  analog  of  coherence.  When 
the  objects  are  not  in  contact  with  each  other,  the  resulting  algorithm  for  a  pairwise 
test  essentially  runs  in  0(1  +  /(M)),  where  f(M)  is  a  function  of  relative  motion  M 


(relative  homogeneous  transformation  matrix  to  describe  the  relative  motion)  which 
the  objects  undergo. 

We  also  analyze  the  complexity  of  our  approaches  to  curved,  models.  The  com- 
plexity of  global  methods  is  0(M3),  where  M  corresponds  to  the  algebraic  complexity 
of  the  resulting  system.  It  is  a  function  of  the  degrees  of  the  polynomials  used  to 
represent  curved  models.  At  each  iteration  the  complexity  of  the  local  methods  is  a 
function  of  the  algebraic  degrees  of  the  boundary  surfaces  and  the  motion  undergone 
by  the  object  pairs. 

The  rest  of  the  paper  is  organized  as  follows.  We  analyze  the  problem  of  interfer- 
ence detection,  review  literature  on  object  models  and  solutions  of  algebraic  systems 
in  Section  2.  In  Section  3,  we  review  the  expected  constant  time  algorithm  for  tracking 
closest  features  between  convex  polytopes  and  show  how  it  can  be  extended  to  inter- 
ference detection  between  general  polyhedral  models.  The  problem  of  closest  features 
and  interference  between  curved  models  are  analyzed  in  Section  4.  In  Section  5,  we 
present  algorithms  between  curved  models  using  equation  solving  approaches.  Finally, 
we  discuss  our  implementation  and  performance  in  Section  6. 

2      Background 

2.1      Coherence  for  Detection  Problem 

The  simplest  algorithms  for  interference  detection  are  based  on  bounding  volumes. 
Typically  they  correspond  to  bounding  boxes,  spheres  or  octrees  etc.  In  case,  two 
objects  are  far  apart,  the  bounding  volumes  approach  can  be  used  to  determine  that 
there  is  no  interference.  These  bounding  boxes  overlap  only  if  the  given  pair  of  objects 
collide  or  are  close  to  each  other.  As  a  result,  any  good  algorithm  for  interference  de- 
tection needs  to  compute  the  spatial  relationship  between  a  pair  of  objects,  when  they 
are  close  to  each  other  and  perhaps  moving  towards  each  other.  Since  no  assumptions 
are  made  on  the  motion  of  the  objects,  it  is  not  possible  to  exactly  compute  the  time  at 
which  the  given  pair  of  objects  would  collide.  Therefore,  in  most  applications  the  col- 
lision detection  routines  are  being  called  frequently.  Due  to  the  nature  of  application, 
there  exists  spatial  and  temporal  coherence  between  successive  instances.  Different 
algorithms  in  the  literature  have  utilized  the  coherence.  Baraff  uses  cached  edges  and 
faces  to  find  a  separating  plane  between  two  convex  polytopes  [3].  However,  this  may 
not  work  when  the  closest  features  between  the  objects  change  or  when  there  is  a  large 
abrupt  motion.  In  [27],  an  algorithm  for  keeping  track  of  closest  features  between  con- 
vex polytopes  using  Voronoi  regions  of  the  features  and  local  geometry  information 
is  used  to  determine  the  minimum  separation  between  them,  thus  to  detect  possible 


collision.  This  works  very  well  for  convex  polytopes.  However,  it  is  difficult  to  apply 
this  approach  directly  to  non-convex  polygonized  curved  models  as  there  may  be  more 
than  pair  of  closest  features  and  the  closest  pair  features  may  correspond  to  curves 
on  the  surface.  In  Section  5,  we  make  use  of  the  coherence  along  with  local  numerical 
methods  and  global  techniques  for  solving  algebraic  systems  to  check  whether  there  is 
a  contact  or  interference  between  such  curved  models. 

2.2  Object  Models 

Polyhedral  models  are  represented  using  boundary  representation  [23].  A  polyhedron 
consists  of  vertices,  edges  and  faces  represented  in  terms  of  parameters  needed  to  de- 
scribe them.  In  addition,  we  also  have  the  topological  information  such  as  adjacencies 
and  incidences  of  all  features  (vertices,  edges,  faces).  Concave  polytopes  are  decom- 
posed into  convex  polytopes  and  represented  by  hierarchical  sub-part  trees.  Curved 
models  consist  of  NURBS  surfaces  and  piecewise  algebraic  surfaces.  NURBS  models 
are  represented  in  terms  of  control  polytopes,  knot  vectors  and  order  continuity  [18]. 
The  convex  hull  of  the  control  polytopes  serves  as  a  bounding  volume  approximation 
for  these  surfaces.  They  are  further  decomposed  into  a  series  of  Bezier  surfaces  using 
knot  insertion  algorithms.  Each  Bezier  surface  corresponds  to  a  rational  parametric 
surface,  which  can  be  represented  in  homogeneous  coordinates  as: 

F(M)  =  (X(s,t),Y(s,t),Z(s,t),W(s,t)).  (1) 

The  Bezier  surface  is  defined  over  the  domain  (s,t)  €  [0, 1]  x  [0, 1].  Algebraic  surfaces 
are  represented  implicitly  as  zero  sets  of  the  form: 

f(x,y,z)  =  o. 

Many  commonly  used  models  like  the  quadrics  or  torus  are  described  using  such 
formulations. 

2.3  Solutions  of  Algebraic  Systems 

It  is  shown  in  Section  5  that  the  problem  of  interference  detection  corresponds  to 
finding  common  solutions  of  a  system  of  algebraic  equations  in  a  particular  domain. 
The  system  may  be  a  zero  dimensional  system  with  a  finite  number  of  solutions  or 
an  overconstrained  system  with  no  common  solutions.  There  are  different  well-known 
methods  for  solving  system  of  polynomials  equations.  Most  of  the  algorithms  in  liter- 
ature are  for  a  system  of  n  equations  in  n  unknowns  represented  as: 

F1(xi,x2,...,xn)     =    0 


F2{xux2,...,xn)    =    0  (2) 

Fn{xi,x2,. .  .,xn)    =    0 

Let  their  degrees  be  rfi,  c?25  •  •  -5  ^n,  respectively. 

Purely  symbolic  approaches  like  Grobner  bases  and  resultants  eliminate  variables 
and  reduce  the  problem  to  finding  roots  of  a  univariate  polynomial.  These  can  be  slow 
even  for  low  degree  systems  and  need  exact  arithmetic.  On  the  other  hand  iterative 
methods  like  the  Gauss-Newton's  method  are  good  for  local  analysis  and  need  good 
initial  guesses  to  all  the  solutions.  In  the  context  of  floating  point  arithmetic,  the  two 
main  approaches  are  those  of  homotopy  methods  [19]  and  matrix  computations  [30]. 
The  latter  algorithms  make  use  of  matrix  formulation  of  resultants  of  polynomial  equa- 
tions and  reduce  the  problem  to  computing  eigenvalues  and  eigenvectors  of  generalized 
companion  matrices.  They  have  been  applied  to  a  number  of  geometric  applications 
and  perform  very  well  in  practice  [30].  The  complexity  of  the  algorithm  is  0(M3), 
where  M  corresponds  to  the  algebraic  complexity  of  the  system.  For  dense  polyno- 
mial systems,  it  is  the  Bezout  bound  corresponding  to  the  product  of  the  degrees  of 
the  equation;  and  for  sparse  polynomial  systems,  it  is  the  BKK  bound  corresponding 
to  the  mixed  volume  of  the  Newton  poly  tope  corresponding  to  each  equations  [6,  41]. 
Homotopy  methods  use  the  solutions  of  a  known  system  along  with  tracing  algorithms 
to  find  the  solutions  of  the  given  system.  The  tracing  steps  corresponds  to  Newton's 
iteration. 

3      Collision  Detection  for  Polyhedra 

In  this  section,  we  summarize  a  simple  and  efficient  collision  detection  algorithm  for 
convex  polyhedra  by  tracking  the  closest  points  between  them  [27].  We  also  extend  it 
to  handle  penetration  between  objects.  The  method  works  by  finding  and  maintaining 
a  pair  of  closest  features  (vertex,  edge,  or  face)  on  the  two  convex  polyhedra,  in  order 
to  calculate  the  Euclidean  distance  between  them  to  detect  possible  collision.  The 
method  is  applicable  in  static  environment,  but  is  especially  well  suited  to  dynamic 
domains  as  the  objects  move  in  a  sequence  of  small,  discrete  steps.  We  take  advantage 
of  the  empirical  fact  that  the  closest  features  change  only  infrequently  as  the  objects 
move  along  finely  discretized  paths.  By  preprocessing  the  polyhedra,  the  algorithm 
runs  in  expected  constant  time  if  the  objects  are  not  moving  very  swiftly.  Even  when 
the  closest  feature  pair  is  changing  rapidly,  the  algorithm  takes  only  slightly  longer 
(runtime  is  proportional  to  the  number  of  feature  pairs  traversed  which  is  a  fuction  of 
the  relative  motion  objects  undergo). 


3.1  Preliminaries 

We  represent  each  convex  polyhedron  using  modified  boundary  representation.  Each 
polyhedron  data  structure  has  fields  for  its  features  (faces,  edges,  and  vertices)  and 
Voronoi  regions  (defined  below). 

Each  feature  (a  vertex,  an  edge,  or  a  face)  is  described  by  its  geometric  param- 
eters and  its  neighboring  features,  i.e.  the  topological  information  of  incidences  and 
adjacencies.  In  addition,  we  also  precompute  the  Voronoi  region  (defined  below)  for 
each  feature  as  well. 

Definition:  A  Voronoi  region  associated  with  each  feature  is  a  set  of  points  closer 
to  that  feature  than  to  any  others.  [38] 

The  Voronoi  regions  form  a  partition  of  space  outside  the  polyhedron  according 
to  the  closest  feature.  The  collection  of  Voronoi  regions  of  each  polyhedron  is  the 
generalized  Voronoi  diagram  of  the  polyhedron.  Note  that  the  generalized  Voronoi 
diagram  of  a  convex  polyhedron  has  linear  size  and  consists  of  polyhedral  regions. 
A  cell  is  the  data  structure  for  a  Voronoi  region.  It  has  a  set  of  constraint  planes 
which  bound  the  Voronoi  region  with  pointers  to  the  neighboring  cells  (which  share 
a  constraint  plane  with  it)  in  its  data  structure.  If  a  point  lies  on  a  constraint  plane, 
then  it  is  equi-distant  from  the  two  features  which  share  this  constraint  plane  in  their 
Voronoi  regions.  Using  the  geometric  properties  of  convex  sets,  applicability  criteria 
(explained  in  Sec. 3. 3)  are  established  based  upon  the  Voronoi  regions.  If  a  point  P 
on  object  A  lies  inside  the  Voronoi  region  of  the  feature  Jb  on  object  5,  then  /g  is  a 
closest  feature  to  the  point  P. 

3.2  Overview 

The  distance  computation  algorithm  which  tracks  the  closest  feature  pair  is  more  fully 
described  in  our  earlier  work  [27]  or  [28].  Here  we  only  give  a  general  overview  of  the 
algorithm. 

Our  method  is  straightforward  in  its  conception.  We  start  with  a  candidate  pair  of 
features,  one  from  each  polyhedron,  and  check  whether  the  closest  points  lie  on  these 
features.  Since  the  objects  (thereby  the  faces  as  well)  are  convex,  this  is  a  local  test 
involving  only  the  neighboring  features  of  the  current  candidate  features.  If  either 
feature  fails  the  test,  we  step  to  a  neighboring  feature  of  one  or  both  candidates,  and 
try  again.  With  some  simple  preprocessing,  we  can  guarantee  that  every  feature  has 
a  constant  number  of  neighboring  features.  This  is  how  we  can  verify  or  update  the 
closest  feature  pair  in  constant  time. 

When  a  pair  of  features  fails  the  test,  the  new  pair  we  choose  is  guaranteed  to  be 
closer  than  the  old  one.    So  when  the  objects  move  and  one  of  the  closest  features 


changes,  we  usually  find  it  after  one  or  two  iterations.  Even  if  the  closest  features  are 
changing  rapidly,  say  once  per  step  along  the  path,  our  algorithm  takes  only  slightly 
longer.  In  this  case,  the  running  time  is  proportional  to  the  number  of  feature  pairs 
traversed  in  this  process.  It  is  not  more  than  the  product  of  the  numbers  of  features  of 
the  two  polyhedra,  because  the  Euclidean  distance  between  feature  pairs  must  always 
decrease  when  a  switch  is  made  [27],  which  makes  cycling  impossible. 

3.3      Applicability  Test 

Given  the  overview  of  the  algorithm,  we  now  present  a  concise  description  of  the  main 
component  of  the  algorithm  —  the  applicability  test  —  which  is  established  based 
upon  Voronoi  region. 

A  Voronoi  region,  associated  with  each  feature,  is  bounded  by  the  constraint  planes 
derived  from  the  feature's  neighbors.  The  applicability  test  is  a  simple  checking  process 
that  verifies  whether  a  point  lies  inside  the  Voronoi  region  of  a  given  feature.  With 
the  preprocessing  procedures  to  guarantee  that  every  feature  has  a  constant  size  of 
neighboring  features,  each  applicability  test  is  only  a  local  test  which  runs  in  0(1)  time. 
When  a  point  fails  the  applicability  test  of  a  given  feature,  the  pointer  associated  with 
each  constraint  plane  provides  a  new,  closer  feature  which  shares  the  same  constraint 
plane  with  the  previous  cell. 

For  a  given  pair  of  features,  feat  a  and  feats-,  on  objects  A  and  B,  we  first  find 
a  pair  of  nearest  points,  Pa  and  Ps,  between  these  two  features.  Then,  we  need  to 
verify  that  feats  is  truly  the  closest  feature  on  B  to  Pa  (i.e.  Pa  lies  inside  the  Voronoi 
region  of  feats)  and  feat  a  is  truly  the  closest  feature  on  A  to  Ps  (i.e.  Ps  satisfies 
the  applicability  test  of  feat  a)-  If  either  check  fails,  a  new  and  closer  feature  (which 
is  usually  one  of  neighboring  features  of  the  previous  features)  is  substituted,  and 
the  new  pair  of  features  is  checked.  Eventually,  we  must  terminate  with  the  closest 
pair,  since  the  distance  between  the  feature  pair  is  strictly  decreasing  through  each 
iteration. 

This  verifying  process  is  demonstrated  in  Fig.l  where  the  previous  closest  feature 
candidate  pair  is  (Fa,Vj,),  Pa  and  VJ,  are  the  two  nearest  points  on  these  two  given 
features.  Though  Pa  satisfies  the  applicability  of  VJ,  (i.e.  Pa  lies  inside  of  H's  Voronoi 
region),  Vj,  fails  Fa's  applicability  test  imposed  by  the  constraint  plane  CP.  Therefore, 
a  new  candidate  pair  (Ea,Vb),  which  is  closer  in  distance,  is  returned.  Then,  the 
algorithm  is  called  again  upon  (Ea,Vb)  to  verify  whether  this  is  the  closest  feature- 
pair  iteratively,  until  the  nearest  points  on  two  features  both  satisfy  the  applicability 
tests  of  each  other.  Then,  the  algorithm  stops  and  returns  the  closest  feature-pair. 
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Figure  1:  Applicability  Test:  7?i  and  R2  are  the  Voronoi  regions  of  Fa  and  £?„  re- 
spectively. (Fa,Vj,)  — >  (^a,Vb)  since  Vj,  fails  the  applicability  test  imposed  by  the 
constraint  plane  CP. 


Figure  2:  Preprocessing  of  a  vertex's  conical  applicability  region  and  a  face's  cylindrical 
applicability  region 


3.4      Subdivision  Procedure 

In  general,  each  face  of  a  typical  convex  polyhedron  has  four  or  five  edges  in  its 
boundary  and  each  vertex  has  three  or  four  incident  edges  .  When  a  face  has  more 
than  5  edges  in  its  boundary  or  a  vertex  has  more  than  5  incident  edges,  its  Voronoi 
region  is  preprocessed  by  subdividing  the  whole  region  into  smaller  cells.  In  Fig. 2,  Ra, 
the  Voronoi  region  of  Va  which  has  8  incident  edges,  is  subdivided  into  jRi,  i?2,  and 
i?3;  Rb,  the  Voronoi  region  of  F\>  which  has  8  edges  in  its  boundary,  is  subdivided  into 
R4,  i?5,  and  Rq.  After  subdivision,  each  of  the  new  cells  has  only  3  or  4  constraint 
planes.  This  subdivision  procedure  is  a  simple  linear  time  routine  which  can  be  done 
as  a  precomputation  step.  It  guarantees  that  when  the  algorithm  starts,  each  Voronoi 
cell  has  a  constant  number  of  constraint  planes.  Consequently,  each  applicability 
test  described  above  runs  in  constant  time  for  updating  the  closest  feature  pairs  in  a 
dynamic  environment. 


3.5  Implementation  Issues 

In  order  to  minimize  online  computation  time,  we  do  all  the  subdivision  procedures 
and  build  all  the  Voronoi  regions  first  as  one-time  computational  cost.  We  do  not  need 
to  reconstruct  all  the  Voronoi  regions  for  each  polyhedron  as  the  objects  move.  We 
only  need  to  transform  the  candidate  features  and  recompute  the  nearest  points,  since 
local  applicability  constraints  are  all  we  need  for  tracking  the  closest  pair  of  features. 
That  is,  we  transform  the  point  coordinates  using  relative  homogeneous  transformation 
matrices  (available  from  either  dynamic  calculations  or  motion  transformations  as 
described  in  Sec  2.2)  and  leave  constraint  plane  equations  fixed.  For  example,  given 
two  objects  moving  in  space,  their  motions  with  respect  to  the  origin  of  the  world  frame 
can  be  characterized  by  the  transformation  matrices  Ta  and  Tg  respectively.  Then, 
their  relative  motion  can  be  represented  by  the  homogeneous  relative  transformation 

Ti  rp—\rp 

AB  —  *B    1A- 

This  algorithm  outputs  the  pair  of  closest  features  and  also  the  distance  between 
two  objects.  If  the  distance  is  less  than  or  equal  to  S  (a  small  safety  margin  defined 
by  the  user  or  determined  by  a  given  environment),  then  the  objects  collide. 

3.6  Hierarchical  Representation  for  Non-Convex  Objects 

We  extend  the  collision  detection  algorithm  for  convex  polyhedra  to  non-convex  ob- 
jects using  hierarchical  representation.  We  assume  that  each  nonconvex  object  is  given 
as  a  union  of  convex  pieces  oris  composed  of  several  nonconvex  subparts,  each  of  these 
can  be  further  represented  as  a  union  of  convex  subparts  or  a  union  of  non-convex 
pieces.  We  use  a  sub-part  hierarchy  tree  to  represent  each  nonconvex  object  (includ- 
ing curved  objects  which  will  be  discussed  later  Fig. 6).  At  each  node  of  the  tree, 
we  store  either  a  convex  sub-part  or  the  union  of  several  convex  subparts.  First,  we 
construct  the  convex  hull  for  each  node  and  work  up  the  tree  as  part  of  preprocessing 
computation.  We  also  include  the  convex  hull  of  the  union  of  sub-parts  in  the  data 
structure.  The  convex  hull  of  each  node  is  the  convex  hull  of  the  union  of  its  children. 
For  instance,  The  root  of  this  sub-part  hierarchy  tree  is  the  nonconvex  object  with  its 
convex  hull  in  its  data  structure. 

At  each  time  step,  we  examine  the  possible  interference  by  a  recursive  algorithm. 
The  algorithm  first  checks  for  collision  between  two  parent  convex  hulls.  Of  course,  if 
there  is  no  interference  between  two  parents,  there  is  no  collision  and  the  algorithm 
stops  and  returns  the  closest  feature  pair  between  two  convex  hulls  of  the  objects.  If 
there  is  a  collision,  then  it  expands  their  children.  If  there  is  also  a  collision  among  the 
leaves,  then  the  algorithm  recursively  proceeds  down  the  tree  to  determine  if  a  collision 
actually  occurs.  In  this  recursive  manner,  the  algorithm  only  signals  a  collision  if  there 
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is  actually  an  impact  among  sub-parts  of  two  objects;  otherwise,  there  is  no  collision 
between  the  two  objects. 

For  complex  objects,  using  a  deep  hierarchy  tree  with  lower  branching  factor  will 
keep  down  the  number  of  nodes  which  need  to  be  expanded.  This  approach  guarantees 
that  we  find  the  earliest  collision  between  two  non-convex  objects. 

3.7     Penetration  Detection  for  Convex  Polyhedra 

The  core  of  the  collision  detection  algorithm  is  built  upon  the  concepts  of  Voronoi 
regions  for  convex  polyhedra.  As  mentioned  earlier  in  Sec  3.2,  the  Voronoi  regions  form 
a  partition  of  space  outside  the  polyhedron.  Therefore,  the  algorithm  can  possibly  run 
into  a  cyclic  loop  when  interpenetration  occurs,  if  no  special  care  is  taken  to  prohibit 
such  events.  Hence,  if  the  polyhedra  can  overlap,  it  is  important  that  we  add  a  module 
which  detects  interpenetration  when  it  occurs  as  well. 

•  Pseudo  Voronoi  Regions: 

This  module  can  be  constructed  based  upon  similar  ideas  of  space  partitioning 
to  the  interior  space  of  the  convex  polyhedra.  The  partitioning  does  not  have  to 
form  the  exact  Voronoi  regions  since  we  are  not  interested  in  knowing  the  closest 
features  between  two  interpenetrating  polyhedra  but  only  to  detect  such  a  case.  A  close 
approximation  with  simple  calculation  can  provide  the  necessary  tools  for  detecting 
overlapping. 

This  is  done  by  partitioning  the  interior  of  a  polyhedron.  We  first  calculate  the 
centroid  of  each  convex  polyhedron,  which  is  the  weighted  average  of  all  the  vertices, 
and  then  construct  the  constraint  planes  of  each  face  to  the  centroid  of  the  polyhedron. 
These  interior  constraint  planes  of  a  face  F  are  the  hyperplanes  passing  through  the 
centroid  and  each  edge  in  F's  boundary  and  the  hyperplane  containing  the  face  F. 
If  all  the  faces  are  equi-distant  from  the  centroid,  these  hyperplanes  form  the  exact 
Voronoi  diagrams  for  the  interior  of  the  polyhedron.  Otherwise,  they  will  provide  a 
reasonable  space  partition  for  the  purpose  of  detecting  interpenetration. 

The  data  structure  of  these  interior  pseudo  Voronoi  regions  is  very  much  like  the 
exterior  Voronoi  regions  described  in  Sec  3.1.  Each  region  associated  with  each  face 
has  e-f  1  hyperplanes  defining  it,  where  e  is  the  number  of  edges  in  the  face's  boundary. 
Each  hyperplane  has  a  pointer  directing  to  a  neighboring  region  where  the  algorithm 
will  step  to  next,  if  the  constraint  imposed  by  this  hyperplane  is  violated.  In  addition, 
a  type  field  is  added  in  the  data  structure  of  a  Voronoi  cell  to  indicate  whether  it  is 
an  interior  (pseudo)  or  exterior  Voronoi  region. 

The  exact  Voronoi  regions  for  the  interior  space  of  the  polyhedron  can  also  be 
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constructed  by  computing  all  the  equi-distant  bisectors  of  all  facets.  But,  such  an 
elaborate  precomputation  is  not  necessary  unless  we  are  also  interested  in  comput- 
ing the  degree  of  interpenetration  for  constructing  some  type  of  penalty  functions  in 
collision  response  or  motion  planning. 

As  two  convex  polyhedral  objects  move  and  pass  through  each  other,  the  change 
of  feature  pairs  also  indicate  such  a  motion  since  the  pointers  associated  with  each 
constraint  plane  keep  track  of  the  pseudo  closest  features  between  two  objects  during 
penetration  phase  as  well.  However,  to  ensure  convergence,  we  will  need  to  construct 
the  exact  Voronoi  regions  of  the  interior  space  or  use  special  case  analysis  to  guaran- 
tee that  each  switch  of  feature  necessarily  decreases  the  distance  between  candidate 
features. 

4     Interference  Detection  for  Curved  Surfaces 

In  this  section,  we  analyze  the  problem  of  interference  detection  between  curved  ob- 
jects represented  as  spline  models  or  piecewise  algebraic  surfaces.  We  show  that  these 
problems  reduce  to  finding  solutions  of  a  system  of  algebraic  equations.  In  particular, 
we  present  algebraic  formulations  corresponding  to  closest  points  determination  and 
geometric  contacts. 

4.1      Closest  Features  Formulation 

Given  the  homogeneous  representation  of  two  parametric  surfaces,  F(s,tf)  =  (X(s,t), 
Y(s,t),  Z(s,t),  W{s,t))  and  G(u,v)  =  (X{u,v),  Y(u,v),Z(u,v),  W(u,v)),  the  clos- 
est features  correspond  to  points  or  curves  on  the  surface.  The  closest  features  are 
characterized  by  the  property  that  the  corresponding  surface  normals  are  collinear. 
This  can  be  expressed  in  terms  of  the  following  variables.  Let 

Fii(s,<,M,u,ai)  =  (F(s,t)  -  G(u,v)) 

Fi2(M,u,u,ai)  =  cti(Gu(u,v)  x  Gv(u,v)) 

F2i(s,£, u,v,a2)  =  (Fs(s,t)  xFt(s,t)) 

F22{s,t,u,v,a2)  =  a2(Gu(u,v)  x  Gv(«,i))), 

where  Fs,Ft,Gu,G„  correspond  to  the  partial  derivatives.  The  closest  features  be- 
tween the  two  surfaces  satisfy  the  following  equations: 

/0\ 
Fi(s,t,u,v,aua2)  =  F1i(s,t,u,v,a1)  -  Fi2(s,t,u,i7,ai)  = 


0 

.0; 


(3) 
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F2(s,  t,  w,  v,  au  a2)  =  F2i(s,  t,  u,  v,  a2)  -  F22(s, i,  u,  u,  a2)  = 


0 


This  results  in  6  equations  in  6  unknowns.  However,  we  are  only  interested  in  the 
solutions  in  the  domain  of  interest  (since  each  surface  is  defined  on  a  subset  of  the 
real  plane).  These  constraints  between  the  closest  features  can  also  be  expressed  as: 

Hi(5,*,u,u)  =(F(M)-G(«,t7))»Gtt(«,t;)=  0 

H2(s,t,u,v)  =(F(s,t)-G(u,v))*Gv{u,v)=  0  (4) 

H3(s,t,w,v)  =  (F(s,t)-G(u,v))*Fa(s,t)  =  0 

H4{s,t,u,v)  =  (F(s,t)  -  G(u,v))  •Ft(s,t)  =  0 

where  •  corresponds  to  the  dot  (or  inner)  product.  This  results  in  four  equations  in 
four  unknowns. 

Let  us  analyze  the  algebraic  complexity  of  these  two  systems  of  equations  cor- 
responding to  closest  features.  Lets  consider  the  first  system  corresponding  to  (3). 
Given  2  rational  parametric  surfaces  F(s,t)  and  G(u,v),  both  of  their  numerators 
and  denominators  are  polynomials  of  degree  n.  The  numerator  and  denominator  of 
Fn(s,2,  u,  v,ai)  have  degree  2n  and  2n  due  to  subtraction  of  two  rational  polynomials. 
As  for  Fi2(s,£,  u,  v,  oji),  the  degrees  of  numerators  and  denominators  of  the  partials 
are  In  —  1  and  2n  respectively  in  the  given  variables  (due  to  quotient  rule).  Taking 
the  cross  product  doubles  the  degrees  for  both  the  numerator  and  denominator;  there- 
fore, the  degrees  for  the  numerator  and  denominator  of  Fi2(s,£,u,  v,a\)  are  An  —  2 
and  4n  respectively.  To  eliminate  a\  from  Fi(s,£,u,i>,ai),  we  get  frfd*'/"'^1))  — 
v1  ,v\s,t,,u'v,Q  W  =  ir/7  )s't,u'v,Q*\  i  After  cross  multiplication  and  clearing  out  the  de- 
nominators,  we  get  two  polynomials  of  degree  12n  —  2  each.  Once  again,  by  the  same 
reasoning  as  stated  above  for  subtraction  and  dot  product  of  rational  polynomials, 
both  the  numerators  and  denominators  of  F2i (5,  <,w,  ^,0-2)  and  F22(s,tf,u,  v,a2)  have 
degrees  of  An  —  2  and  An.  By  similar  method  mentioned  above,  we  can  eliminate  a2 
from  F2(s,i,  w,  v,a2).  We  get  two  polynomial  equations  of  degree  16n  —  4  each  after 
cross  multiplication.  As  a  result  the  system  has  a  Bezout  bound  of  (12n  — 2)2(16n— 4)2. 

Each  equation  in  the  second  system  of  equations  has  degree  An  —  1  (obtained  after 
computing  the  partials  and  subtraction  of  two  rational  polynomials)  and  therefore  the 
overall  algebraic  complexity  corresponding  to  the  Bezout  bound  is  (An  —  l)4.  Since 
the  later  system  4  results  in  a  lower  degree  bound,  in  the  rest  of  the  analysis  we  will 
use  this  system.  However,  we  are  only  interested  in  the  solutions  in  the  domain  of 
interest  (since  each  surface  is  defined  on  a  subset  of  the  real  plane). 
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Figure  3:  Tangential  intersection  and  boundary  intersection  between  two  Bezier  sur- 
faces 


Baraff  has  used  the  formulation  in  Eqn.  3  to  keep  track  of  closest  points  between 
closed  convex  surfaces  [3]  based  on  local  optimization  routines.  The  main  problem  is 
finding  a  solution  to  these  equations  (3)  for  the  initial  configuration.  In  general,  these 
equations  can  have  more  than  one  solution  in  the  associated  domain  (even  though  there 
is  only  one  closest  point  pair)  and  the  optimization  routine  may  not  converge  to  the 
right  solution.  A  simple  example  is  the  formulation  for  the  problem  of  interference 
detection  between  two  spheres.  There  is  only  one  pair  of  closest  points,  however 
Equations  (4)  have  four  pairs  of  real  solutions. 

Given  two  algebraic  surfaces,  f(x,y,z)  =  0  and  g(x,y,z)  =  0,  the  problem  of 
closest  features  determination  can  be  reduced  to  finding  roots  of  the  following  system 
of  6  algebraic  equations: 


/(zi,3/i,zi) 

#(*2,t/2,22) 
fy{Xl,yi,Zl) 

^  fz{x\,y\,zx)  j 


=    0 
=    0 

=    «1 


(   9x(x2,y2,Z2)    \ 

gy{x2,y2,Z2) 
^  gz(x2,y2,z2)  ) 


(5) 
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(  x1  \ 


fx2\ 


=    a2 


(  gx(x2,y2,z2)  ^ 

gy{x2,y2,z2) 
\  9z{x2,y2,z2)  j 


Given  two  algebraic  surfaces  of  degree  n,  we  can  eliminate  a!  by  setting  ■MXl'Vl'*1)  = 
Jv[xi,yi,zi)  _  fz[xi,yi,zi)     ^fter  cross  multiplication,  we  have  a  polynomial  equation  of 

9y  (*2 .1/2  ,Z2)  5z(*2, J/2,22)  r  '  f        J  n 

2n  —  2,  since  each  partial  has  degree  of  n  —  1  and  the  multiplication  results  in  the  degree 
sum  of  2n  —  2.   Similarly,  to  eliminate  a2,  we  set      fl ~X2    \  =      ,yi~y2    ,  =       *i~*2 

J '  "  Si(^2.y2,«2)  5»(^2,y2,Z2)  5z(x2.J/2.22) 

and  the  degree  of  the  resulting  polynomial  equations  is  n.  We  have  six  quations 
after  eliminating  ai  and  a2:  two  of  degrees  (2n  —  2)  and  four  of  degrees  n  respec- 
tively (2  from  eliminating  a2  and  2  from  /(xi,t/i,2i)  and  g(x2,y2,z2)).  Therefore, 
the  Bezout  bound  of  the  resulting  system  can  be  as  high  as  iV  =  (2n  —  2)2n4.  In 
general,  if  the  system  of  equations  is  sparse,  we  can  get  a  tight  bound  with  Bern- 
stein bound  [5].  The  Bernstein  bound  for  Eqn.  5  is  n2(n2  +  3)(n  —  l)2.  Canny 
and  Emiris  calculate  the  Bernstein  bounds  by  using  sparse  mixed  resultant  formula- 
tion [6].  For  example,  the  Bernstein  bounds1  for  the  case  of  n  =  2,3,4,5,6,7,8,9 
are  28,432,2736,11200,35100,91728,210112,435456,  while  the  Bezout  bounds  are 
64, 1296,9216,40000, 129600,345744,  •  •  •  respectively.  Even  for  small  values  of  n,  the 
bound  on  the  number  of  solutions  can  be  large,  and  therefore  the  algebraic  complexity 
of  computing  the  closest  points  can  be  extremely  high,  even  using  the  exact  tight 
bound.  In  our  applications  we  are  only  interested  in  the  real  solutions  to  these  equa- 
tions in  the  corresponding  domain  of  interest.  The  actual  number  of  real  solutions 
may  change  as  the  two  objects  undergo  motion  and  some  configurations  can  result  in 
infinite  solutions  (e.g.  when  a  closest  pair  corresponds  to  a  curve  on  each  surface,  as 
shown  for  two  cylinders  in  Fig.  4.  ).  As  a  result,  it  is  fairly  non-trivial  to  keep  track 
of  all  the  closest  features  between  objects  and  updating  them  as  the  objects  undergo 
motion. 

4.2      Contact  Formulation 

The  problem  of  interference  detection  corresponds  to  determining  whether  there  is 
any  contact  between  the  two  objects.  In  particular,  it  is  assumed  that  in  the  be- 
ginning two  objects  are  not  overlapping.  As  they  undergo  motion,  we  are  interested 
in  knowing  whether  there  is  any  precise  contact  between  the  objects.  There  are  two 
types  of  contacts.  They  are  tangential  contact  and  boundary  contact.  In  this  section, 
we  formulate  both  of  these  contact  determination  problems  in  terms  of  a  system  of 


1These  figures  are  calculated  by  John  Canny  and  Ioannis  Emiris  using  their  code  based  on  the 
sparse  mixed  resultant  formulation  [6]. 
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algebraic  equations.  In  the  next  section,  we  describe  how  the  algorithm  tests  for  these 
conditions  as  the  object  undergoes  rigid  motion. 

•  Tangential  Intersection  :  This  corresponds  to  a  tangential  intersection  between 
the  two  surfaces  at  a  geometric  contact  point,  as  in  Fig. 3(a).  The  contact  point 
lies  on  the  interior  of  each  surface  (as  opposed  to  being  on  the  boundary  curve) 
and  the  normal  vectors  at  that  point  are  collinear.  These  constraints  can  be 
formulated  as: 

F(s,t)  =  G(u,v) 
(Fa{s,t)xFt(s,t))*Gu(u,v)  =  0  '    (6) 

(F.(M)  xFt(s,t))*Gv(u,v)  =  0 

The  first  vector  equation  corresponds  to  a  contact  between  the  two  surfaces  and 
the  last  two  equations  represent  the  fact  that  their  normals  are  collinear.  They 
are  expressed  as  scalar  triple  product  of  the  vectors.  The  last  vector  equation 
represented  in  terms  of  cross  product  corresponds  to  three  scalar  equations. 
We  obtain  5  equations  in  4  unknowns.  This  is  an  overconstrained  system  and 
has  a  solution  only  when  the  two  surfaces  are  touching  each  other  tangentially. 
However,  we  solve  the  problem  by  computing  all  the  solutions  to  the  first  four 
equations  using  global  methods  and  substitute  them  into  the  fifth  equation.  If 
the  given  equations  have  a  common  solution,  than  one  of  the  solution  of  the  first 
four  equation  will  satisfy  the  fifth  equation  as  well.  For  the  first  three  equations, 
after  cross  multiplication  we  get  3  polynomial  equations  of  degree  2n  each.  The 
dot  product  results  in  the  addition  of  degrees  of  the  numerator  polynomials  (each 
of  these  partials  has  numerator  polynomial  of  degree  (2n  —  1)).  Therefore,  we 
get  a  polynomial  equation  of  degree  6n  —  3  from  the  fourth  equation.  Therefore, 
the  Bezout  bound  of  the  system  corresponding  to  the  first  four  equations  is 
bounded  by  N  =  (2n)3(6n  — 3),  where  n  is  the  parametric  degree  of  each  surface. 
Similarly  for  two  algebraic  surfaces,  the  problem  of  tangential  intersection  can 
be  formulated  as: 


1  fx(x,y,z)  ^ 
fv(x,y,z) 

\  fz(x,y,z)  ) 


=  a 


f(x,y,z)  =  0 

g(x,y,z)  =  0  (7) 

'  9x(x,y,z)  N 

9y(x,y,z) 
\  9z(x,y,z)  j 


In  this  case,  we  obtain  4  equations  in  3  unknowns  (after  eliminating  a)  and 
these  equations  correspond  to  an  overconstrained  system  as  well.    These  over- 
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Figure  4:  Closest  features  between  two  different  orientations  of  a  cylinder 

constrained  system  is  solved  in  a  manner  similar  to  that  of  parametric  surfaces. 
The  Bezout  bound  for  the  first  three  equations  is  M  =  n2(2n  —  2). 

•  Boundary  Intersection  :  Such  intersections  lie  on  the  boundary  curve  of  one  of 
the  two  surfaces.  Say  we  are  given  a  Bezier  surface,  defined  over  the  domain, 
(s,2)  G  [0, 1]  x  [0, 1],  we  obtain  the  boundary  curves  by  substituting  s  or  t  to  be 
0  or  1.  The  resulting  problem  reduces  to  solving  the  equations: 


F(j,l)  =  G(u,tO 


(8) 


Other  possible  boundary  intersections  can  be  computed  in  a  similar  manner.  The 
intersection  points  can  be  easily  computed  using  global  methods.  An  example 
has  been  shown  in  Figure  3(b) 

Two  objects  have  a  contact  if  one  of  these  sets  of  equations,  ((6)  or  (8))  for  para- 
metric surfaces  and  (7)  for  algebraic  surfaces,  have  a  common  solution  in  their  domain. 

In  a  few  degenerate  cases,  it  is  possible  that  the  system  of  equations  (6)  and  (8) 
have  an  infinite  number  of  solutions.  One  such  example  is  shown  for  two  cylinders  in 
Fig. 4.  In  this  case  the  geometric  contact  corresponds  to  a  curve  on  each  surface,  as 
opposed  to  a  point.  These  cases  can  be  detected  using  resultant  methods  as  well  [30]. 

4.3      Penetration  Analysis 

In  most  examples,  it  is  difficult  to  express  the  motion  as  a  closed  form  function  of 
time.  As  a  result,  it  is  almost  impossible  for  the  algorithm  to  compute  the  exact  time 
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to  collision  as  at  instance  we  only  have  the  information  about  the  object's  position 
and  geometry.  Based  on  the  algorithm  for  polyhedral  models  and  curved  models 
(presented  in  the  next  section)  it  is  possible  for  us  to  know  that  the  two  objects  are 
getting  closer.  As  a  result,  the  main  routine  of  the  system  can  check  for  interference 
detection  at  shorter  time  steps.  However,  it  is  possible  that  at  a  given  time  instance 
the  two  objects  are  not  colliding  and  at  the  next  instance  they  have  penetrated  into 
each  other.  The  algebraic  constraints  for  contact  analysis,  based  on  tangential  contact 
or  boundary  contact,  are  satisfied  only  if  the  two  objects  have  a  tangential  or  boundary 
contact.  They  may  no  longer  be  satisfied,  when  the  two  objects  have  penetrated.  In 
this  section,  we  formulate  conditions  for  penetration  between  curved  models  and  how 
they  can  be  used  for  collision  detection. 

Two  objects  are  penetrating,  if  there  is  a  finite  area  of  contact  between  them 
(as  shown  in  Fig.  5).  In  contact  analysis,  we  have  assumed  that  there  is  a  point 
contact  (or  a  higher  order  contact  along  the  boundary  or  on  the  tangent  plane).  The 
penetration  therefore  results  in  a  one-dimensional  intersection  curve  corresponding  to 
the  intersection  of  the  two  surface  boundaries.  In  this  section,  we  show  that  there  is 
still  a  real  solution  to  the  subset  of  equations  corresponding  to  (6),  (8)  or  (7). 

Lets  consider  boundary  intersections  as  in  Fig.  5(b).  Lets  assume  that  there  was  a 
contact  along  the  boundary  curve  corresponding  to  F(s,  1)  and  at  the  given  instance 
the  two  surfaces  have  penetrated.  As  a  result,  if  we  compute  the  common  solutions  of 
the  equations: 

F(*,l)  =  G(ti,tO  (9) 

there  is  a  real  solution  in  the  corresponding  domain,  which  may  not  correspond  to 
the  boundary  curve  of  G(u,v)  or  there  are  more  than  one  distinct  solutions  to  these 
equations.  However,  any  small  perturbation  along  the  boundary  curve  F(s,l)  would 
result  in  a  solution  to  these  equations  in  the  corresponding  domain  and  based  on  that 
we  can  detect  interference. 

Now,  lets  consider  the  tangential  contact  between  two  algebraic  surfaces  as  in 
Fig.  5(a).  Similar  analysis  would  apply  to  the  tangential  intersection  of  piecewise 
parametric  surfaces  as  well.  To  check  for  tangential  contact,  we  use  the  equation 
solving  methods  to  compute  the  solutions  of  the  following  four  equations  in  four  un- 
knowns: 

f{x,y,z)  =  0 

g(x,y,z)  =  Q  (10) 

/  fx(x,y,z)  \_a(  9x{x,y,z)  \ 

\   fy{x,y,z)   )  V  9y(x,V,z)   J  ' 


18 


Tangent    Plane 


Boundary  Intersection 


overlap  curves 


(a) 


(b) 


Figure  5:  (a)  penetration  resulted  from  tangential  intersection;  (b)  overlap  resulted 
from  boundary  intersection 

The  solutions  are  substituted  into  the  fifth  equation,  as  shown  in  (6)  to  check  for 
contact  analysis.  The  penetration  information  is  obtained  based  on  following  theorem. 
Lets  assume  that  the  intersection  between  the  two  surfaces  is  non-planar. 

Theorem  4.1  Two  algebraic  surfaces,  f(x,y,z)  =  0  and  g(x,y,z)  =  0  have  a  non- 
planar  interpenetration,  iff  there  is  at  least  one  real  solution  of  the  system  of  equations 
(10). 

Proof:  Lets  take  the  case  when  the  two  objects  have  no  geometric  contact  between 
them.  As  a  result,  there  is  no  common  real  solution  to  the  first  two  equations, 
/(x,  y,z)  =  0  and  g(x,  y,  z)  =  0. 

Lets  consider  the  case,  when  the  objects  penetrate.  Moreover,  the  intersection 
is  non-planar.  As  a  result,  the  two  surfaces  intersect  along  a  one  dimensional  space 
curve.  In  a  few  degenerate  cases,  the  intersection  curve  is  higher  dimensional.  Lets 
denote  the  real  curve  of  intersection  as  C .  C  is  an  algebraic  space  curve  corresponding 
to  all  the  common  solutions  of  f(x,y,z)  =  0  and  g(x,y,z)  =  0.  Typically,  C  consists 
of  one  closed  component  in  space  (as  shown  in  Fig. 5).  In  some  cases,  corresponding 
to  non-convex  models  it  may  consist  of  multiple  components.  Eliminating  a  from  the 
last  two  equations  of  (10)  results  in 

G{x,  y,  z)  =  /x(ar,  y,  z)gy(x,  y,  z)  -  fy(x,  y,  z)gx(x,  y,  z)  =  0. 
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The  normal  at  a  point  on  the  intersection  curve  is  obtained  by  taking  the  cross 

product  of  the  two  surfaces  at  that  point.    In  other  points,  the  normal  function  is 

defined  as: 

f  fx{x,y,z)  \       I  gx(x,y,z)  ^ 

N(x,y,z)=       fy(x,y,z)       x       gy(x,y,z)      . 

\  fz(x,y,z)  )       \  gz(x,y,z)  j 

It  follows  from  this  definition  that  G(x,y,z)  corresponds  to  the  Z-component  of  the 
normal  vector  of  the  intersection  curve.  As  a  result,  the  real  solutions  of  the  system 
of  equation  (10)  correspond  to  all  those  points  on  the  intersection  curve  where  the  Z- 
component  of  the  normal  vector  is  zero.  Since  C  consists  of  non-planar  closed  loops, 
there  are  at  least  two  such  points  on  each  component  of  C.  As  a  result,  the  given 
system  of  equations  have  a  real  solution.  Q.E.D. 

The  penetration  criterion  highlighted  in  Theorem  4.1  assumes  that  the  intersection 
curve  is  non-planar.  This  is  typically  the  case  for  most  models  (unless  one  of  them  is 
a  plane  or  if  the  two  models  correspond  to  quadric  surface).  In  fact,  the  non-planarity 
condition  for  penetration  can  be  relaxed  and  it  is  required  that  the  curve  should 
have  a  point  where  the  Z-component  of  the  normal  vector  is  zero.  This  condition 
is  not  satisfied,  if  the  intersection  curve  lies  in  a  plane  parallel  to  the  XY-p\a,ne. 
To  circumvent  this  case,  we  take  three  generic  linear  combinations  of  the  last  three 
equations  highlighted  in  (7)  and  pick  two  of  the  resulting  three  equations  along  with 
f(x,y,z)  =  0  and  g(x,y,z)  =  0  to  formulate  a  new  system.  For  almost  all  linear 
combinations  the  penetration  condition  highlighted  in  Theorem  4.1  is  exact.  Similar 
analysis  can  be  applied  to  the  tangential  intersection  of  two  parametric  surfaces.  In 
case  the  two  models  penetrate  and  there  is  no  solution  for  the  equations  corresponding 
to  boundary  intersection  in  the  given  domain,  there  is  a  real  solution  to  the  first  five 
equations  highlighted  in  (6)  in  the  given  domain.  Again  this  may  only  fail  if  the 
intersection  curve  is  planar  and  there  is  no  point  where  the  Z-component  of  the  normal 
vector  vanishes.  To  improve  the  robustness  of  the  approach  we  replace  the  last  three 
equations  in  (6)  by  their  generic  linear  combinations. 

5      Coherence  for  interference  detection  between 
curved  objects 

In  most  dynamic  environments,  the  closest  features  between  two  moving  objects 
change  infrequently  between  two  time  frames.  We  have  used  this  coherence  property, 
utilizing  spatial  as  well  as  temporal  coherence,  in  designing  the  expected  constant 
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Figure  6:  Hierarchical  representation  of  a  torus  composed  of  Bezier  surfaces 

time  algorithm  for  collision  detection  among  convex  polytopes  and  extended  to  non- 
convex  polytopes  as  well.  In  this  section  we  utilize  this  coherence  property  along  with 
the  algebraic  formulations  presented  in  the  previous  section  for  contact  determination 
between  curved  models. 

Given  two  B-spline  models,  we  decompose  them  into  a  series  of  Bezier  surfaces 
using  knot  insertion  algorithm  [18].  After  decomposition  we  use  a  hierarchical  repre- 
sentation for  the  curved  models.  The  height  of  the  resulting  tree  is  two.  Each  leaf  node 
corresponds  to  the  Bezier  surface.  The  nodes  at  the  first  level  of  the  tree  correspond 
to  the  convex  hull  of  the  control  polytope  of  each  Bezier  surface.  The  root  of  the  tree 
represents  the  union  of  the  convex  hull  of  all  these  polytopes.  The  torus  is  composed 
of  biquadratic  Bezier  surfaces,  shown  at  the  leaf  nodes  (Fig.  6).  The  convex  polytopes 
at  the  first  level  are  the  convex  hull  of  control  polytopes  of  each  Bezier  surface  and 
the  root  is  a  convex  hull  of  the  union  of  all  control  points. 
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The  algorithm  proceeds  on  the  hierarchical  representation,  as  explained  in  Section 
3.  However,  each  leaf  node  is  a  Bezier  surface.  The  fact  that  each  surface  is  non-convex 
implies  that  the  closest  features  of  the  polytopes  may  not  be  a  good  approximation  to 
the  closest  points  on  the  surface.  Moreover,  two  such  surfaces  can  have  more  than  one 
closest  feature  at  any  time.  Given  two  algebraic  surfaces,  we  initially  enclose  them  in 
a  bounding  box  and  test  the  bounding  boxes  for  collision.  At  the  leaf  nodes,  we  check 
the  two  models  for  precise  contact,  once  their  control  polytopes  collide. 

The  problem  of  collision  detection  between  two  Bezier  surfaces  is  solved  by  finding 
all  the  solutions  to  the  equations  (6)  and  (8).  If  the  two  models  are  described  implicitly, 
the  first  collision  corresponds  to  the  solution  of  (7).  A  real  solution  in  the  domain 
to  those  equations  implies  a  geometric  collision  and  a  precise  contact  between  the 
models.  The  algebraic  method  based  on  resultants  and  eigenvalues  is  used  to  find  all 
the  solutions  to  the  equations  (6), (8)  or  (7)  [30].  This  global  root  finder  is  used  when 
the  control  polytopes  of  two  Bezier  surfaces  collide.  At  that  instant  the  two  surfaces 
may  or  may  not  have  a  geometric  contact.  It  is  possible  that  all  the  solutions  to 
these  equations  are  complex.  The  set  of  equations  in  (6)  represents  an  overconstraint 
system  and  may  have  no  solution  in  the  complex  domain  as  well.  However,  we  apply 
the  algebraic  method  to  the  first  four  equations  in  (6)  and  compute  all  the  solutions. 

The  total  number  of  solutions  of  a  system  of  equations  is  bounded  by  the  Bezout 
bound.  The  resultant  method  computes  all  these  solutions.  As  the  objects  move, 
we  update  the  coefficients  of  these  equations  based  on  the  rigid  motion.  We  obtain 
a  new  set  of  equations  corresponding  to  (6)  and  (8),  whose  coefficients  are  slightly 
different  as  compared  to  the  previous  set.  All  the  roots  of  the  new  set  of  equations 
are  updated  using  Gauss-Newton's  method.  The  previous  set  of  roots  are  used  as 
initial  guesses.  Gauss-Newton's  method  works  well,  if  the  relative  motion  between  the 
two  objects  is  small.  This  approach  is  similar  to  homotopy  methods  for  tracing  paths 
to  compute  solutions  of  a  given  system  of  polynomial  equations.  In  case  the  Gauss- 
Newton's  method  does  not  converge  to  a  solution,  we  may  again  use  the  problem 
in  terms  of  an  eigenvalue  problem  and  use  the  previous  eigenvalues  as  a  guess  and 
use  inverse  power  iteration.  This  approach  has  been  explained  in  detail  in  [30].  The 
convergence  of  inverse  power  iteration  is  well  understood  and  it  converges  to  the 
eigenvalue  closest  to  the  guess.  As  a  result,  we  are  performing  local  computations  at 
each  step  corresponding  to  a  few  Gauss-Newton's  iterations  or  inverse  power  iterations. 
This  procedure  represents  an  algebraic  analog  of  the  geometric  coherence  exploited  in 
the  earlier  section. 

Typically,  if  the  objects  have  a  tangential  contact  and  as  the  two  objects  are  moving 
closer  to  each  other,  the  imaginary  components  of  some  of  the  roots  start  decreasing. 
This  involves  the  use  of  complex  arithmetic  while  updating  the  solution  with  Gauss- 
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Newton's  method  or  inverse  power  iterations.  Finally,  a  real  collision  occurs,  when  the 
imaginary  component  of  one  of  the  roots  becomes  zero.  A  real  solution  of  the  resulting 
system  implies  a  collision.  In  case  it  satisfies  all  the  equations  simultaneously,  there 
is  a  tangential  contact  otherwise  there  is  penetration.  The  objects  have  a  boundary 
intersection,  if  the  equations  corresponding  to  boundary  intersection  have  a  real  solu- 
tion, which  is  not  in  the  domain.  As  the  objects  undergo  motion,  the  solution  moves 
and  based  on  the  roots  of  the  given  system  of  equations,  we  can  determine  whether 
there  is  a  boundary  contact  or  penetration. 

It  is  possible  for  the  objects  to  have  multiple  contacts.  In  such  cases  two  distinct 
sets  of  solution  paths  converge  to  real  solutions  in  the  corresponding  domain.  As  a 
result,  the  algorithm  can  easily  find  all  the  multiple  point  of  contacts. 

The  performance  of  the  algorithm  is  dependent  on  the  complexity  of  the  global 
root  finder  and  the  number  of  paths  we  need  to  trace.  The  total  number  of  paths 
correspond  to  the  Bezout  bound  of  the  given  system  and  can  be  fairly  high  for  a  given 
system  of  equations.  Tracing  each  path  at  each  instance  can  be  time  consuming. 
The  choice  of  paths  is  made  based  on  the  solutions  of  the  global  equation  solver. 
In  most  applications,  the  global  root  finder  is  used  when  the  two  objects  are  fairly 
close  to  each  other.  This  happens  when  the  bounding  boxes  corresponding  to  control 
poly  topes  overlap.  As  a  result,  we  only  choose  those  solutions  as  the  start  paths 
whose  imaginary  components  are  relatively  small  or  whose  real  solutions  are  close  to 
the  corresponding  domain.  This  is  a  function  of  the  degree  of  the  resulting  surfaces. 
For  higher  degree  models,  we  pick  initial  solutions  with  high  imaginary  components 
as  well. 

Example  5.1  Lets  illustrate  the  algorithm  on  two  tori  undergoing  motion  and  even- 
tually colliding  at  a  point.  We  assume  that  each  torus  is  enclosed  by  a  bounding  box. 
Such  a  bounding  box  can  either  be  obtained  using  representation  of  the  torus  as  a  union 
of  rational  Bezier  surfaces  (taking  the  convex  hull  of  associated  control  points)  or  from 
the  geometric  definition.  In  particular,  a  torus  can  be  described  as  a  union  of  points 
equi-distant  from  a  circle  in  space.  As  a  result,  a  torus  can  be  described  in  terms  of 
a  circle  C  in  3  space  and  a  radius  r.  The  bounding  box  is  computed  by  enclosing  the 
circle  C  with  the  square  of  minimum  area,  say  S ,  and  its  height  is  r  along  each  side. 
A  canonical  representation  of  the  torus,  assuming  a  circle  of  radius  k  centered  at  the 
origin  is 

F(x,y,z)  =  (y/^  +  y2  -  k)2  +  z2  -  r  =  0. 

In  Figure  7,  we  have  shown  two  tori  moving  towards  each  other.  The  initial  position 
of  the  first  torus  is  described  using  equation 

F(x,  y,  z)  =  64  -  20x2  +  x4  -  20y2  +  2x2y2  +  y4  +  I6z2  +  2x2z2  +  2y2z2  +  z4 
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and  the  initial  configuration  of  the  second  torus  corresponds  to 

G(x,  y,  z)  =  392.0625  -  291x  +  84.5x2  -  12x3  +  x4  +  32.5y2  -  12xy2  +  2x2y2  + 
y4  -  221. bz  +  84x2  -  14x2*  -  14y2z  +  81.5a:2  -  12xz2  +  2x2z2  +  2y2z2  -  14z3  +  z4 


To  check  for  contact  we  add  the  following  three  equations  corresponding  to  the  first 
order  differentials  of  F(x,y,z)  and  G(x,y,z): 


-40x  +  4x3  +  4xy2  +  4xz2 

-40y  +  4x2y  +  4V3  +  4yz2 

32z  +  4x2  z  +  4y2  z  +  4z3 


—  a 


-291  +  169x  -  36x2  +  4x3  -  12y2  +  4xy2  +  84z  -  28. xz  -  12z2  +  4xz2 

65y  -  24xy  +  4x2y  +  4J/3  -  28yz  +  4yz2 
-227.5  +  84x  -  14x2  -  14y2  +  163z  -  24xz  +  4x2z  +  4y2z  -  42z2  +  4z3 


To  solve  for  the  equations,  we  take  linear  combinations  of  these  three  equations, 
eliminate  a  and  reduce  the  problem  of  contact  determination  to  a  system  of  4  equations 
in  3  unknowns.  As  the  objects  undergo  motion,  these  system  of  equations  is  updated 
in  the  following  manner.  Let  the  rigid  motion  of  the  first  object  be  represented  as 


+  Tlt 


where  R\  is  an  orthogonal  matrix  corresponding  to  the  rotation  and  T\  corresponds  to 
the  translation.   This  can  be  expressed  as 


(x\ 

(  x\ 

y 

=  Ri 

y 

\J  ) 

\  z  ) 

(  x) 

(x\ 

y 

=  J?( 

y 

\  z  J 

\J  I 

-Ti) 


The  corresponding  position  of  the  new  object  is  obtained  by  substituting  for  (x,y,z) 
in  F(x,y,z)  and  obtaining  the  equation  F(x,y,z)  =  0  corresponding  to  the  position 
of  the  torus  after  undergoing  motion.  Similarly,  we  compute  a  new  representation 
G(x,y,~z)  ==  0  for  the  second  torus.  Notice  that,  it  undergoes  a  motion  represented 
by  i?2  and  T2.  Given  these  equations,  we  formulate  the  problem  of  contact  determi- 
nation by  taking  their  differentials  and  taking  the  linear  combinations  of  the  resulting 
equations.  The  solution  of  the  previous  set  of  equations  are  used  as  the  starting  guesses. 
This  process  is  repeated  until  there  is  a  contact  or  penetration  between  the  two  objects 
or  their  bounding  boxes  do  not  overlap  anymore. 
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6  Implementation  and  Performance 

We  have  implemented  the  algorithms  described  in  this  paper  and  tested  their  per- 
formance on  many  cases.  The  collision  detection  routine  for  convex  polytopes  is  very 
efficient  and  gives  us  a  constant  time  performance  of  50  to  70  usee  (x  10-6  sec)  per  pair 
of  objects  (independent  of  object  complexity)  on  SGI  Indigo2  XZ.  The  hierarchical 
approach  works  well  for  concave  objects,  given  a  good  convex  decomposition  of  the 
object.  We  are  able  to  perform  collision  detection  in  real  time  using  our  polyhedral 
collision  detection  algorithm. 

Using  the  spatial  and  temporal  coherence,  we  can  keep  track  of  the  closest  feature 
pairs  between  convex  curved  objects  and  the  control  polytopes  of  non-convex  curved 
objects  efficiently  as  well.  In  Fig. 7(a)  -  7(f)  (from  left  to  right,  top  to  bottom)  we 
demonstrate  the  algorithm  on  two  tori.  Each  torus  is  represented  hierarchically  as 
shown  in  Fig. 6.  In  particular,  Fig. 7(c)  and  Fig. 7(d)  represent  the  same  time  instant. 
In  Fig. 7(c)  the  convex  hulls  of  two  tori  collide  but  in  Fig. 7(d)  the  control  polytopes 
of  the  subparts  do  not  touch  each  other.  Once  two  control  polytopes  of  the  subparts 
also  collide  Fig.7(e),  we  use  the  global  methods  to  find  the  initial  solutions  of  (6) 
and  (8)  and  followed  by  Gauss-Newton  method  to  update  the  solutions  at  each  time 
instance.  It  involves  complex  arithmetic.  Although  the  hierarchical  representation  is 
in  terms  of  biquadratic  Bezier  patches,  we  use  the  implicit  representation  of  a  torus 
(as  that  results  in  a  system  with  lower  Bezout  bound).  Each  torus  is  represented  as 
an  algebraic  surface  of  degree  4.  Its  equation  is  obtained  by  squaring  the  expressions 
in 

F(x,  y,  z)  =  (yW*/2  -  k)2  +  z2  -  r  =  0. 

Since  the  torus  is  a  closed  object,  there  are  no  boundary  intersections  and  we  keep  track 
of  tangential  intersections  for  geometric  contacts.  The  resulting  system  of  equations 
corresponding  to  (6)  are  of  degrees  4,4  and  6.  The  Bezout  bound  of  this  system  is 
96.  Using  resultants,  the  problem  is  reduced  to  finding  the  eigen-decomposition  of  a 
96  x  96  matrix  [30].  The  eigen-decomposition  takes  slightly  more  than  a  second  on 
the  IBM  RS/6000  and  all  the  solutions  at  that  instant  are  complex.  As  the  tori  move 
towards  each  other,  we  eventually  track  8  of  these  solutions,  whose  imaginary  parts 
are  decreasing.  This  uses  Gauss-Newton's  method  at  each  time  instance.  Finally  a 
real  solution  is  obtained  corresponding  to  the  instance  shown  in  Fig. 7(f). 

7  Conclusion 

In  this  paper  we  have  described  efficient  algorithms  for  contact  determination  in  dy- 
namic environments.    The  algorithms  are  based  on  spatial  and  temporal  coherence 
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between  successive  instances  and  applied  to  polyhedral  models  and  curved  objects 
whose  boundaries  can  be  described  as  piecewise  algebraic  sets.  This  include  NURBS 
surfaces  commonly  used  in  computer  graphics  and  geometric  modeling.  The  resulting 
algorithm  make  use  of  hierarchical  representations,  along  with  local  numerical  and 
global  algebraic  methods  for  equation  solving.  These  algorithms  have  been  imple- 
mented and  work  well  in  practice. 

Many  issues  related  to  contact  determination  for  general  dynamic  environments 
still  remain  open.  In  an  environment  with  N  objects,  we  want  to  avoid  the  0(N2) 
pairwise  tests  at  each  instance.  Some  techniques  to  speed  up  the  run  time  based 
on  uniform  spatial  subdivision  are  proposed  in  [42]  and  algorithms  of  complexity 
0(Nlog2N)  have  been  highlighted  for  spheres  in  [25]  and  axis-aligned  bounding  boxes 
in  [16].  However,  they  are  suitable  for  static  environments  but  do  not  take  advantage 
of  coherence  in  dynamic  environments.  Furthermore,  there  are  no  good  algorithms  for 
contact  determination  between  deformable  models. 
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Figure  7:  (a)  -  (d)  Collision  detection  algorithm  applied  to  two  tori 
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Figure  7:  (e)  -  (f)  Collision  detection  algorithm  applied  to  two  tori 
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